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ABSTRACT 

We present results of fully three-dimensional magnetohydrodynamic (MHD) simulations of disk accre- 
tion to a slowly rotating magnetized star with its dipole moment inclined at an angle O to the rotation 
axis of the disk (which is assumed to be aligned with the spin axis of the star). The main goal was to inves- 
tigate the pattern of magnetospheric flow and the disk-star interaction for a variety of inclination angles 
O. We observed that at O = 0°, the disk stops at magnetospheric radius r^, and matter flows to the star 
through axisymmetric funnel flows, as observed in earlier axisymmetric simulations. However, when the 
dipole moment of the star is inclined, then the flow becomes non- axisymmetric. The non-axisymmetry 
becomes notable at very small inclination angles O ~ 2° — 5°. The pattern of magnetospheric flow is 
different at different 8. For relatively small angles, O < 30°, the densest matter flows to the star mostly 
in two streams, which follow paths to the closest magnetic pole. The streams typically co-rotate with 
the star, but they may precess about the star for 8 <; 10°. At intermediate angles, 30° ^ 8 < 60°, the 
streams become more complicated, and often split into four streams. For even larger angles, 8 ^ 60°, 
matter accretes in two streams, but their geometry is different from the streams at small 8. Magnetic 
braking changes the structure of the inner regions of the disk. It creates a region of lower density (a 
"gap") for <^ 4rm. A ring of higher density forms at r ~ r,„ for 6 <^ 30°. For r <^ [2 — 3)rm, the 

azimuthal velocities are sub-Keplerian. The inner region of the disk at r ~ is warped. The warping 
is due to the tendency of matter to co-rotate with inclined magnetosphere. The normal of the inner 
warped part of the disk is close to the magnetic axis of the dipole. The accreting matter brings positive 
angular momentum to the (slowly rotating) star tending to spin it up. The corresponding torque Nz 
depends only weakly on 8. The angular momentum flux to the star near the star's surface is transported 
predominantly by the magnetic field; the matter component contributes ~ 1% of the total flux. The 
torques and Ny are also calculated and these may give a slow precession of the symmetry axis of the 
star. The angle 8 was fixed in simulations because the time scale of its evolution is much longer that 
that of the simulations. Results of simulations are important for understanding the nature of classical T 
Tauri stars, cataclysmic variables, and X-ray pulsars. These stars often show complicated spectral and 
photometric variability patterns, which may be connected with the structure of magnetospheric flows. 
The magnetospheric structure of stars with different 8 can give different variability patterns in observed 
light curves. This can provide information about inclination angles 8 in different stars. A notable result 
of the present simulations is the formation of multiple streams in the accretion flows near the star for 
intermediate inclination angles. This may give short-scale quasi-periodic variability in the light curves 
of some stars. 



L INTRODUCTION 

1.1. Observations. Disk accretion to a rotating star with 
a dipole magnetic field is important in a number of objects, 
including classical T Tauri stars, cataclysmic variables and 
X-ray binaries. Commonly, the magnetic axis of the dipole 
is not aligned with the rotation axis of the disk. The disk 



accretion to rotating neutron stars with misaligned dipole 
fields gives rise to X-ray pulsars in binary systems (see, 
e.g., Triimper et al. 1985; Sheffer et al. 1992; Bildsten 
et al. 1997; Deeter et al. 1998). The periodic variability 
in some cataclysmic variables is thought to be due to the 
rotation of a star with a misaligned dipole field (e.g., Livio 
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& Pringlc 1992; Warner 1995; Warner 2000; Wickramas- 
inghe, Wu, & Ferrario 1991). A similar cause may give 
the photometric and spectral variability of some classical 
T Tauri stars (e.g., Herbst et al. 1986; Bouvier & Bertout 
1989; Johns & Basri 1995; Hartmann 1998; Bouvier, et al. 
1999, 2003; Petrov, et al. 2001). 

Disk accretion to an inclined dipole rotator occurs in dif- 
ferent systems. For specificity we adopt a scaling of the pa- 
rameters appropriate for classical T Tauri Stars (CTTSs). 
Application to other objects such as neutron stars may be 
done by rescaling the parameters. 

Following the discovery of variability of classical T Tauri 
stars (CTTSs) (Herbst et al. 1986; Bertout et al. 1988), 
it was suggested that this variability may be due to rota- 
tional modulation of emission of "hot spots" on the sur- 
face of the star which are connected with magnetospheric 
funnel flows (Camenzind 1990; Konigl 1991). Many ob- 
servations of CTTSs support this model. In some cases 
(e.g., AA Tau), the variability may be due to the warp 
of the inner part of the disk (Bouvier et al. 1999, 2003). 
Photometric and spectral variability of CTTSs often shows 
complicated patterns which have not been explained. In 
order to understand the nature of the variability of differ- 
ent CTTSs, it is important to calculate the detailed MHD 
flow pattern in the magnetosphere outside of the inclined 
rotator. The magnetic field of CTTSs may not be a pure 
dipole (e.g., Safier 1998, Smirnov et al. 2003). However, 
this paper considers the case where the star's field is a 
dipole field. 

1.2. Theory of Disk Accretion to an Inclined Rotator. 
The problem of disk accretion to an inclined dipole ro- 
tator has been investigated theoretically by a number of 
authors. These works considered the the possible warping 
of the inner part of the disk. 

It was noticed in 1970s that in case of diamagnetic disk 
(when magnetic field does not penetrate the disk) the mag- 
netic pressure will be different above and below the disk 
so that there is a net force perpendicular to the disk. This 
force causes warping of the disk (Scharlemann 1978; Aly 
1980). This early work considered the case in which the ro- 
tational axis of the distant disk Jl^ and spin of the star fl^ 
coincide. Later, the analysis was extended to the case in 
which the CI4 and $7* are in different directions (Lipunov & 
Shakura 1980; Lipunov, Semenov & Shakura 1981; Horn & 
Kundt 1989; Lai 1999). The net force acting on the inner 
regions of the disk (averaged over the rotation of the star) 
acts so as to make the normal to the warped disk perpen- 
dicular to both the magnetic axis of the dipole and spin of 
the star. The disk is expected to precess around the spin 
axis of the star. Subsequently, Lai (1999) extended this 
analysis to cases of imperfectly conducting disks. Later, it 
was shown that viscosity may decrease warping (Horn & 
Kundt 1989; Lai 1999). Agapitou, Papaloizou & Terquem 
(1997) and Terquem & Papaloizou (2000) conclude that 
bending waves may also lead to warping, and the ampli- 
tude of the warp will depend on viscosity parameter a. 

The theoretical models typically involve a number of 
simplifying approximations. Predictions of these models 
should be checked numerically. This paper, however, does 
not focus on this topic, and no special simulations were 
arranged to enhance the warping. 

No theories were done on the main topic of this paper 



- analysis of magnetospheric flows in the inclined magne- 
tosphere. This problem is principally three-cUmensional 
and needs full 3D simulations. Some of the important 
problems, however, were investigated analytically and nu- 
merically in an axisymmetric approach, which is valuable 
for subsequent 3D analysis. 

2.3. 2D Simulations of the Disk-Star Interaction. Dif- 
ferent aspects of the disk-star interaction for the case of 
an aligned rotator were investigated with two-dimenionsal, 
axisymmtric MHD simulations (Hayashi, Shibata & Mat- 
sumoto 1996; Miller & Stone 1997; Hirose et al. 1997). 
Special attention was given to the problem of opening 
of the coronal magnetic field lines and associated out- 
flows (Goodson, Bohm & Winglee 1999; Fend & Elstner 
2000). This process was earlier investigated theoretically 
by Lovelace, Romanova, & Bisnovatyi-Kogan (1995), and 
Bardou & Heyvaerts (1996). Much less attention has been 
given to numerical investigation of the dipole interaction 
with the inner regions of the disk and modelling of mag- 
netospheric flows, which were investigated analytically by 
Ghosh & Lamb (1979a,b), Shu et al. (1994), Ostriker & 
Shu (1995), Wang (1995), Li & Wickramasinghe (1997), 
Koldoba et al. (2002a), and others. This problem was not 
fully investigated by MHD simulations because of the dif- 
flculty of creating consistent "quiescent initial conditions." 
In most simulations the initial magnetic braking of the disk 
by the magnetospheric field leads to the fast accretion of 
disk matter to the star with velocity close to free-fall (e.g., 
Hayashi et al. 1996; Miller & Stone 1997). Furthermore, 
the strong initial twist forms at the boundary between the 
disk and corona, and this acts to initiate outflows. 

Romanova et al. (2002, hereafter R02) developed con- 
sistent quiescent initial conditions, which lead to slow, vis- 
cous time-scale accretion avoiding the free-fall accretion of 
earlier studies. This opened up the possibility of investi- 
gating the physics of the disk-star interaction and magne- 
tospheric flows in great detail on the viscous time-scale. 
Many aspects of the earlier predicted theories were recon- 
sidered (R02). It was shown that (1) the magnetic field 
lines in the inner regions of the disk have a tendency to be 
closed, while external lines may be open (see also recent 
paper by Kiiker, Henning, Riidiger 2003); (2) the inner 
regions of the disk experience magnetic braking and their 
angular velocity is smaller /larger than Keplerian (depend- 
ing on angular velocity of the star); (3) the magnetic field 
is dominant in the spinning-up/spinning-down of the star; 
(4) the main force lifting matter up to the funnel flow is 
the pressure gradient force, while magnetic force is negli- 
gibly small. Different aspects of earlier proposed models 
were confirmed or rejected. With detailed axisymmetric 
analysis in hand, it was easier to under take the three- 
dimensional MHD simulations of the disk accretion to the 
inclined dipole. 

In order to understand accretion to a misaligned rota- 
tor, we created a special three-dimensional MHD simula- 
tion code based on a "cubed sphere" grid. This grid has 
a number of advantages over spherical or Cartesian grids 
(Koldoba et al. 2002b, hereafter K02). A similar grid was 
developed for the surface of the sphere for geophysical ap- 
phcations (see, e.g., Ronchi, lacono, & Paolucci 1996). In 
contrast with these authors, we perform simulations in the 
three-space and use as a base the Godunov-type numer- 
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Fig. 1. — Initial conditions for simulations. A star of radius ii* is threaded with the dipole magnetic field with the 
magnetic momentum fi inclined relative to the z— axis at an angle 6. i?max - is the size of the simulations region. The 
solid lines with arrows represent magnetic field lines. The level lines of the disk are shown. The inner radius of the disk 
is located at R^. 



ical scheme which is similar to one described by Powell 

ct al. (1999). Wc performed sample simulations of three- 
dimensional accretion to an inclined rotator for an incli- 
nation angle O = 30° and obtained the first interesting 
results. We observed that after one rotation of the inner 
radius of the disk, the initial axisymmetry of the disk is 
destroyed, and matter flows to the star along two streams. 

In this paper we present the results for the full range of 
inclination angles G. The main goals of the simulations are 
the following: (1) to investigate the structure of the mag- 
netospheric flows at different inclination angles 8, which 
is important for subsequent analysis of spectral and photo- 
metric variability of CTTSs and other stars; (2) to inves- 
tigate the properties of the disk surrounding the magneto- 
sphere; (3) to investigate the rate of spinning up/down of 
the star and its dependence on 0; and (4) to investigate 
how inclination angles G may change with time. In future 
papers we will show detailed analysis of the 3D streams 
and hot spots, and will derive the light curves from simu- 
lations (Romanova et al. 2004a). 

In §2, the numerical model is described. In §3, test sim- 
ulations are presented. In §4, the main numerical results 
for different G are discussed. In §5, the magnetic brak- 
ing of the disk is analyzed. In §6, the fluxes of matter 
and torque to the star arc calculated. In §7, the main 
conclusions of this work are presented and the possible 
observational consequences are discussed. 



2. NUMERICAL MODEL 

We investigate the case in which the angular velocity of 
the star $7 = $7, and that of the disk are ahgned and in 
the z direction. The magnetic moment fj, is inclined rela- 
tive to z axis by an angle G as shown in Figure 1. We fix 
inclination angle G and do not change it during simula- 
tions, thus reflecting the fact that the angular momentum 
of the star is much larger than the angular momentum 
of accreting matter. In general, the angular momenta of 
the star and that of the disk may be misaligned, but in 
the present paper we consider only the aligned case. The 
initial magnetic field is considered to be a dipole field, 
B = 3(/x • R)R/i?^ — /x/i?^. We use a reference frame 
(X, Y, Z) rotating with the star with the Z axis aligned 
with the star's rotation axis. This rotating frame is orien- 
tated such that /Lx is in the {X, Z) plane. 

We solve the system of ideal MHD equations. We de- 
compose the magnetic field B into the "main" dipole com- 
ponent of the star Bq, and the component, Bi, induced 
by currents in the disk and in the corona. This helps 
to minimize errors in calculating the magnetic force (e.g., 
Tanaka 1994). In the case of the inclined rotating dipole 
the dipole moment changes with time. It rotates with an- 
gular velocity CI so that the "main" field Bq also changes 
with time. Consequently in the induction equation there 
is a large term involving Bq. To overcome this difficulty 
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we use a coordinate system rotating with angular veloc- 
ity r2, in which the magnetic moment of the star fi and 
the "main" field Bq do not depend on time. The system 
of MHD equations in the coordinate system rotating with 
the star is 

^^+V-(pv)=0, (1) 



dt 



dt 

+ V • T = pg + 2p V X n - p 12 X (n X R) , (2) 



d{ps) 



V • (psv) = , 



dt 

-=Vx(vxB) 



(3) 
(4) 



(K02) where p is the plasma density, s is the entropy, Tik 
is the stress tensor with Tif. — p6ik + pViVk + (B^(5ifc/2 — 
BiBk)/ATT. Here v is the velocity of plasma in the rotat- 
ing frame. The transformation to the inertial (observer's) 
frame is u = v + n x R. One can see that equation (2), 
which is the Eulcr equation in the rotating frame, has new 
terms compared to the usual Euler equation, namely, the 
Coriolis and centrifugal forces. The Godunov-type numer- 
ical code was used in the "cubed sphere" coordinate sys- 
tem. The code was tested and verified to give the correct 
solution for Bondi accretion to a gravitating object with a 
monopolc magnetic field (K02). 

2.1. Reference Values and Example for T Tauri Stars 

The reference value for distance is denoted Rq, and it 
is taken to be the initial value of the inner radius of the 
disk. Thus the dimensionless inner disk radius is iJ^ = 1 
initially. We take the dimensionless radius of the star to 
be R^ = 0.35. The reference value for the velocity is taken 
to be the Keplerian velocity at Rq, vq = {GM/RoY^"^ . 
The reference angular rotation rate is ujq = vq/Rq, and 
the corresponding time-scale is to = Ra/vo- We also use 
the rotation period at r = i?o, Pa = 27rto. For a given 
magnetic field strength at Rq we can define a reference 
density as pq = /vq and reference pressure pq = pv^ . A 
reference magnetic moment of the star is then po = BqRq- 
Thus, the calculated variables are: R' = R/Rq, t' = t/to, 
p' = p/po, v' = v/vo, B' = B/Bo, p' = p/po- We also in- 
troduce the dimensionless time, T = t/Po- Subsequently, 
the primes are dropped, but the physical values can easily 
be restored for a particular case. 

Here, we discuss the numerical parameters for a typ- 
ical T Tauri star. We take the mass and radius of the 
star to be M = O.SMq and i?* = 1.8Rq. Our reference 
length is approximately equal to the inner radius of the 
disk Ro « 2.86i?* « 3.6 x lO" cm. The reference velocity 
is Wo ~ 1.93 X 10^ cm/s, and the corresponding time-scale 
is to ~ Ro/vo ~ 1-89 X 10'' s. The size of the simula- 
tion region corresponds to i?max = 0.34 AU. The period 
of Keplerian rotation of the inner radius of the disk is 
Pq « 1.38 days. 

Consider B^, = lO'"* G at the surface of the star. 
Then at R = Rq, the reference magnetic field is Bq = 
B,(i?*/i?o)^ ~ 42.7 G and the reference magnetic moment 
is po = BqRo « 2.0 X 10^'' Gcm~^. The reference density 
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cm 



is po = 4.89 X 10-12 g/cm3 or n = 3.06 x 10 
which is typical for T Tauri star disks. The reference 
mass accretion rate is Mq = pqVqR^ « 1.2 x 10^^ g/s « 
1.92 X 10"'' Mg/year. 



2.2. Initial and Boundary Conditions 

To derive the initial distribution of density and pres- 
sure in the disk and corona, we use quiescent initial con- 
ditions which were developed for our axisymmetric two- 
dimensional (R02) and test three-dimensional simulations 
(K02). The main idea for creation of quiescent initial con- 
ditions was: (1) to minimize the initial magnetic braking 
between the slowly rotating dipole and matter of the disk 
rotating with Keplerian velocity, and (2) to minimize the 
effect of the rapid twisting of the magnetic field between 
the disk and corona as a result of the fast turn on of the 
disk rotation. The first process typically leads to fast ac- 
cretion of matter in the disk with a velocity close to free- 
fall, while the second process leads to the formation of a 
strong B(j) component near the disk and a transient matter 
outburst to the corona (e.g., Hayashi et al. 1996; Miller & 
Stone 1997). To minimize these factors, we had the corona 
rotate with the Keplerian velocity of the disk at some dis- 
tance r > Rd from the axis. Thus, the corona is differen- 
tially rotating with constant angular velocity on cylinders 
of equal radius r, oj{r) = const (see also Romanova et al. 
1998). 

We establish the pressure balance on the boundary 
between the disk and corona such that the disk below 
this boundary has high density pd and low temperature 
Td, while corona above this boundary has low density 
pc « pd and high temperature Tc >> Td- Then we 
derive the density and pressure distribution in the whole 
region taking into account the differential rotation of the 
disk and corona, and force balance between gravitational, 
pressure, and centrifugal forces (see R02 for details). In 
the typical density distribution, density in the disk grad- 
ually increases outward, while density in the corona de- 
creases while moving to the axis (see Figure 1 of R02). 
We take the ratios: p^/pd = 0.03 and Td/Tc = 0.03. 
Note, that in two-dimensional simulations we took the 
ratios, Pc/ pd = 0.01 and Td/Tc = 0.01. We have cho- 
sen softer ratios for three-dimensional simulations, in or- 
der to get a smaller Alfven velocity va = B/^Airp, and 
larger time-step At ~ AR/va- Test two-dimensional and 
three-dimensional simulations with both sets of ratios have 
shown that in both cases the corona does not influence the 
disk flow appreciably. 

We place the initial inner radius of the disk Rd at the 
magnetospheric radius rm , where the magnetic pressure of 
the dipole field balances the matter pressure of the disk, 
p + pv"^ = B^/Stt. We chose Bq at this point and pd such 
that Rd = r^n = This choice is based on the fact that we 
often observe that at small inclination angles the streams 
always start from the magnetospheric radius rm (see also 
Pringle & Rees 1972; Ghosh & Lamb 1979a,b for theoret- 
ical base). This equilibrium was checked numerically in 
two- and three-dimensional simulations (see R02, and §3.1 
of this paper). 

We suppose that the star rotates relatively slowly in the 
sense that its corot ation radiu s is at rcor = 3. This corre- 
sponds to = ^/GM/r^.„. w 0.19 = OMUk*, where 
GM = 1 in our dimensionless units. For the case of 
T Tauri stars, this corresponds to a slow rotation rate 
(T « 9.4 days for the parameters used in §2.1 and in R02). 
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Higher rotation rates including the "propeller" regime will 
be described in Romanova et al. (2004b). 

The boundary conditions arc similar to those used in 
K02. At the inner boundary R = i?*, we take "free" 
boundary conditions d/dR = for all variables. This 
results in the "absorption" of the incoming matter so that 
there is no standoff shock. Further, the inner boundary is 
treated as a perfect conductor rotating at the rate ft = ilz. 
In the reference frame rotating with the star the flow ve- 
locity is parallel to B at i? = i?* . The boundary condition 
at iJ, on the magnetic field has d{RB^)/dR = 0. At the 
outer boundary R = i?max, free boundary conditions are 
taken for all variables. Wc investigated the possible in- 
fluence of the outer boundary conditions by running cases 
where the simulation regions had i?max = 4.9, 14, and 68. 
We find the same results except for the case i?max = 4.9 
where the accretion rate decreases too fast because the 
reservoir of matter in the disk is too small. Results for the 
medium and large regions are very close, so that we take 
.Rmax = 14 as a standard size of the simulation region. 

2.3. The Grid and Parallelization of the Code 

As mentioned the inner boundary of the simulation re- 
gion (the "star") is taken to be at i?* = 0.35. The size of 
external boundary depends on the grid (see below). 

The spherical grid was inhomogeneous in the 
i?— direction. The inhomogeneity was such that cells 
at any distance R were approximately square with 
Ai? ~ R/S.9 for fixed A0. This grid gives a high space 
resolution close to the star which is important in this 
problem. At the same time, simulations in a very large 
region may be performed. The inhomogeneous grid is 
a smooth analog of nested grids used for example by 
GBW99. 

The angular resolution was taken to be Nq = 29 for each 
block in most cases, and Ng = 41 in test cases. The num- 
ber of points in R— direction determines the size of the 
simulation region. The main simulations were done with 
the radial grid Nr = 70 which corresponds to i?max = 
14 PS 40i?H. « 0.34 AU. Simulations at the larger grid 
Nr = 96 corresponding to ii^ax = 68 « 194i?* « 1.6 AU 
were also performed for testing the external boundary con- 
ditions. 

Note, that each of the six blocks consists of 70 x 29 x 29 
cells. In comparison with a spherical coordinate system, 
this grid corresponds to Nr x Ng x = 70 x 58 x 116. 
This grid gives high resolution in the vicinity of the star 
and inner region of the disk. For the typical grid Ng = 29, 
the size of the smallest cells near the surface of the star is 
Ai? = 0.019 w 0.053i?,. This cell size is 1.7 times larger 
than in our two-dimensional simulations of the correspond- 
ing problem (R02), but about 1.8 times smaller than cells 
in the innermost grid of GBW99 in their two-dimensional 
simulations, who had resolution Ai? w O.li?* in the inner- 
most grid. At larger distances from the star, the grid size 
becomes larger but at the magnetospheric distance i? = 1 
it is still very small, Ai? « O.li?*, and continues to be suf- 
ficiently small in the inner region of the disk i? <J 5 which 
is the most important for our current investigation. At 
the external boundary, Rmax = 14, Ai? = 2i?*, which is 
small compared with the disk thickness at this distance. 
These numbers show that this resolution and grid are sat- 



isfactory for three-dimensional simulations of the matter 
flow aroimd the magnetized star. Note, that in K02 it was 
shown that even sparser grids such as 26 x 15 x 15 give 
satisfactory results. 

Fast runs of the three-dimensional code at high grid res- 
olutions become possible due to parallelization of the code. 
First, the code was naturally parallelized to 6 processors, 
where each of six blocks was assigned to one of six proces- 
sors. Next, we divided each block into equidistant layers 
Niayers in dircction. Division onto Niayers = 2, 3, 4, ... 

etc. allowed us to use Nproc = 6 X Niayers — 12, 18, 24, ... 

number of processors. We typically used Niayers = 8, that 
is, 48 processors. Simulations with larger numbers of pro- 
cessors were also possible and still produced an increase 
in the computation speed, but the efficiency of using the 
parallel machine at Nproc > 48 declined. 

The reported simulations were performed mainly on the 
fastest V2 nodes of the "Velocity Cluster" at the Cornell 
Theory Center. Each of the V2 nodes is a dual proces- 
sor 2.4 GHz Pentium 4 computer. Typical simulations of 
T = lOPo rotations with the grid 70 x 29 x 29 for each of 
6 blocks took approximately 40-50 wall-clock hours on the 
48 processors of V2 nodes. This was a reasonable time of 
simulations in that it allowed us to investigate the three- 
dimensional flows at different inclination angles (this pa- 
per) and other parameters (Romanova et al. 2004b), and 
to perform test runs with the grid 70 x 41 x 41. 

3. TEST NUMERICAL SIMULATIONS 

First, we discuss tests of the code with the above men- 
tioned initial and boundary conditions in hydrodynamic 
simulations, that is, B = 0. Secondly, we discuss tests of 
the code for accretion to a rotating star with an aligned 
dipole magnetic field, = 0. Thirdly, we compare re- 
sults from our three-dimensional simulations with analo- 
gous two-dimensional axisymmetric simulations. 

3.1. Hydrodynamic Simulations (B = 0) 

Hydrodynamic simulations are important for checking 
the initial setup of the problem. Simulations were done 
for grids A^i = 70 x 29 x 29 and = 70 x 41 x 41 for each 
of six blocks. The first grid was the base grid for subse- 
quent runs with the magnetic field, while the second grid 
was used for comparisons. Simulations with both grids 
have shown that the disk-star system stays in equilibrium 
during many rotations T > 15Po, and only slightly evolves 
due to numerical viscosity. Inward velocity associated with 
numerical viscosity is very small, v « 0.001 — O.OOSvKd, 
where VKd = \/ GM/Rd is the Keplerian velocity at the 
distance r = R^. For longer times, the outer regions of the 
disk become thicker because the grid resolution in that re- 
gion is relatively low. No high velocities were observed in 
the corona. The initial differential rotation in the corona 
was gradually disturbed, and the corona become almost 
homogeneous but still low-density. The accretion rate to 
the star was calculated to be M = 0.01 — 0.02 in the di- 
mensionless units (see §2.2) for both grids A''i and iV2. 
Note that the accretion rate obtained subsequently in the 
magnetic cases is about 10 times larger. 

Thus, hydrodynamic simulations have shown that the 
disk-star system stays in the gravity/pressure/centrifugal 
equilibrium for a sufficiently long time. In three- 
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Fig. 2. — The figure shows results of three-dimensional MHD simulations of disk accretion to an aligned rotator (0 = 0) 
for T = and T = 5. Time T is measured in units of Keplerian rotation period at the inner radius of the disk r = Rd = 1- 
The isodensity contours (background) are shown at p = 0.15. The solid lines represent the magnetic field lines. The grid 
70 X 29 X 29 is taken for each of six sectors of the cubed sphere grid. Only part of the simulation region i? < 3 is shown. 



dimensional simulations we were able to follow the evo- 
lution during 5—12 rotations Pq. Compared to our two- 
dimensional simulations (R02) which were followed until 
T > 50 rotation, here, we do not have viscosity incorpo- 
rated into the code. Longer simulations with incorporated 
viscosity will be done in the future research. 

In 2D simulations we noticed, that the magnetospheric 
funnel flow establishes during one rotation period of the 
inner disk. Subsequently, the flow changes insignificantly 
during the subsequent 50 rotations and variation depends 
on the incoming matter flux. In simulations of the accre- 
tion to the inclined rotator, we observed that the flow also 
establishes during the first several rotations, so that as 
presented in this paper simulations which lasted 5—12 ro- 
tations arc valuable for understanding of 3D funnel flows. 

3.2. Accretion to an Aligned Rotator, = 

Next, we performed three-dimensional simulations of ac- 
cretion to a magnetized star. First, we chose the case 
where the magnetic axis is perpendicular to the disk and 
G = 0. We took the standard parameters for initial and 
boundary conditions as described in §2.1. We observed 
that during the first two rotation periods, matter of the 
disk moved inward, but later it stopped at the magne- 
tospheric radius « 1 — 1.2 and went up/down out the 
disk plane, forming magnetospheric funnel fiows. The fiow 
is cylindrically-symmetric, as it should be. The shape of 
the funnel flow is similar starting from T w 3 and during 
larger T. Figure 2 shows magnetospheric flow at T = 
and T — b. Here and in subsequent plots we show only 
the inner part of the simulation region, i? < 3, in order 
to resolve the funnel fiow and the inner region of the disk. 
The magnetospheric flow at = has a relatively low 
density, so that we show isocontours at the density level 
p = 0.15. 

The star rotates slowly, 17, = 0.04i7,x with the co- 
rotation radius at rcor = 3. Thus, the magnetic held lines 
threading the disk at r ^ 3 have an angular velocity lower 



than the Keplerian, and they transfer negative angular mo- 
mentum to the disk matter. This leads to the slow inward 
flow of matter to the star with velocity v « 0.03 — 0.05uifd. 
This flow results from the magnetic braking, but the ve- 
locity of the flow is much smaller than that obtained us- 
ing non-equilibrium initial conditions (e.g., Hayashi et al. 
1996). 

To understand the behavior of the magnetic fleld, we 
chose a set of magnetic field lines which initially threaded 
the disk at a distance r 2.8, and show the same set of 
lines at subsequent moments of time. One can see that ac- 
creting matter drags these field lines inward. This demon- 
strates that magnetic field lines are well frozen to the disk. 
They rotate with the disk in the azimuthal direction and 
slowly move inward in the radial direction. No opening 
of magnetic field lines was observed, because the corona 
right above the disk is matter-dominated (e.g., Romanova 
et al. 1998). 

We should note that the boundary between the disk and 
magnetosphere of the aligned dipole is a possible site for 
onset of 3D instabilities, such as the Relegh- Taylor insta- 
bility which may lead to direct accretion of matter through 
magnetosphere to the star (see Arons & Lea 1976a, b, 
Scharlemann 1978). No accretion through 3D instabilities 
was observed in our simulations. Possibly special simula- 
tions are needed with even higher resolution. 

3.3. Comparison of Two-dimensional and 
Three-dimensional Simulations 

We compared results of axisymmetric three-dimensional 
simulations with corresponding two-dimensional simula- 
tions. Two-dimensional simulations were done in spherical 
coordinates with the grid x Ng = 70 x 58 which is the 
same as the cubed sphere grid Nr x x = 70 x 29 x 29 for 
each block. Figure 3 shows the meridional cross-sections of 
matter flow for both cases after T — 5 rotations. Only the 
inner region R < 4.5 is shown in order to resolve the matter 
flow around the magnetosphere. One can see that results 
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Fig. 3. — Comparison of the density profiles for disk accretion to an aligned rotator obtained from two-dimensional, 
axisymmetric simulations in spherical coordinates (R02) (top panel), and from three-dimensional simulations in the 
cubed-sphere coordinates (bottom panel), both at T = 5. 




Fig. 4. — Simulations of accretion to an inclined dipole rotator for an inclination angle — 5°. Background represents 
density at the level p — 0.15. Time T is measured in periods of Keplerian rotation at i? = 1. 



of two-dimensional and three-dimensional simulations are 
qualitatively similar. However, in 3D simulations the mag- 
netosphere radius is ~ 10 — 20% larger than in 2D case: the 
funnel flow starts at « 0.9— 1.1 in the two-dimensional 
simulations and at r,„ w 1.0 — 1.2 in the three-dimensional 
simulations. Besides, the funnel flow in 3D simulations 



forms later, after two rotations, compared to one rotation 
in 2D simulations (R02). These differences seems to be 
connected with the difference in the accretion rate. Accre- 
tion rate in 3D case is about 10 — 20% smaller than in 
2D case: Aho ~ 0.11 - 0.13 versus Nho ~ 0.12 - 0.16. 
In present simulations the inward accretion of matter oc- 




Fig. 5. — Simulations of disk accretion to a misaligned dipole rotator for the case where the dipole moment is inclined 
by an angle 6 = 5°. The background represents the density surface where p — 0.35. Time T is measured in periods of 
Keplerian rotation at i? = 1. 



curs as a result of small initial magnetic braking (see §2.2). 
The magnetic braking may have a different nature in 2D 
and 3D cases which may possibly lead to the observed dif- 
ference between results in 2D and 3D simulations. In 2D 
simulations this braking is possibly stronger, which leads 
to the enhanced accretion compared to 3D simulations. 

Here, we should point out that magnetospheric structure 
does depend on the accretion rate, but does not depend 
on the mechanism of accretion and angular momentum 
transport. In R02 and in recent 3D simulations (R03), we 
introduced the simplified a— type viscosity term to the 
equations, which permitted us to vary the accretion rate. 



In real accretion disks the transport of the angular momen- 
tum is connected most probably with magneto-rotational 
instability (e.g., Balbus & Hawley 1991, 1998). Numerical 
investigation of this instability needs high space resolution 
in the disk and typically only part of the disk is considered 
(e.g., Hawley, Gammie & Balbus 1995; Brandenburg et al. 
1995; Stone, et al. 1996; Fleming, Stone, & Hawley 2000). 
In the present simulations the high grid resolution was set 
in the vicinity of the dipole, while resolution father out in 
the disk was not sufficient to resolve this instability. 

These tests provide assurance about the reliability of 
our three-dimensional MHD code. It is important that at 
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Fig. 6. — Matter flow around an inclinded dipole rotator with the dipole inchnation angle Q 
represents the density isocontours at p = 0.35. 



15°. The background 



8 = 0° the code supports axisymmetry of the flow. We 
next investigate disk accretion to an inclined dipole rota- 
tor. 

4. RESULTS OF SIMULATIONS AT DIFFERENT 
INCLINATION ANGLES 

We performed a number of simulations for different in- 
clination angles, 6 = 2°, 9 = 5°, 6 = 15°, 6 = 30°, 

9 = 45°, e = 60°, and 9 = 75°. We separately discuss 
results for small, medium and large inclination angles. 



4.1. Small Inclination Angles: — 2° 
9 = 15° 
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First, it is important to investigate the difference be- 
tween matter flow at zero, 9 = 0°, and very small incli- 
nation angles. We observed that magnetospheric flow be- 
comes axisymmetric starting from very small inclination 
angles, 9 = 2° and 9 = 5°. At such angles, at low density 
the "windows" of even lower density form in the magneto- 
spheric flow, while at larger densities, matter accretes in 
two streams. 

Figure 4 shows initial state (left panel) and magneto- 
spheric flow at T = 5 (right panel) for relatively low den- 
sity p — 0.15. One can see that compared to the case with 
9 = 0° (Figure 2), the windows of lower density form. 
Windows precess about the rotation axis. The plots are 
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Fig. 7. — The figure shows different projections of the flow for an inclination angle 8 = 15° at time T = 6. The 
background represents the density a,t p — 0.35; the solid lines are magnetic field lines. In panel (a) the fi vector is directed 
towards the observer; in panel (b) the orientation is the same as in Figure 6; panel (c) shows the view along the Z axis 
with vector /i. directed downward; panel (d) shows the side view along the Y axis. 



shown in the coordinate system rotating with the star, and 
the fi vector of the dipole is fixed in the {X, Z) plane and 
has the same direction in each of the two plots of Figure 
4. Precession of the windows means that its rotation is 
different from that of the star. Both the star and the disk 
rotate counterclockwise looking from the +z pole towards 
the disk. 

The windows rotate in the same direction but faster 
which means that the period of precession is shorter than 
that of the star. The faster rotation is due to the fact that 
inner regions of the disk near the magnetosphere rotate 
with an angular velocity close to Keplcrian uj « ^Kir — 
1) « 1 while the star rotates slowly with w 0.19. The 
windows rotate with an angular velocity between these two 
values. 

There is of course a range of densities in the magneto- 
spheric flow. At higher densities, p ^ 0.15, the windows 
become larger, and at even higher densities, p ^ 0.3 — 0.4 
matter flows in dense narrow streams. These streams are 
expected to be important in the occultation of stellar light 
from an observer. They also form hot spots where they 
hit the stellar surface. 

The area at the star covered by the hot spots depends 
on density. At p « 0.15, hot spots cover Agpots ~ 40% of 



the surface of the star, while at p « 0.3 — 0.4, they cover 
Aspots « 20 — 30% of the surface. At larger densities of 
the stream, the area is even smaller (see detailed analysis 
in R03). 

Subsequently, we show the high density stream compo- 
nent, but point out that flow pattern is different at the 
different density levels. 

Figure 5 shows the accretion flow at the density level 
p = 0.35. One can see that after T — 2 rotations two 
streams form. The azimuthal width of the streams is sev- 
eral times larger than their thickness in the poloidal direc- 
tion. The streams precess around the Z axis in the same 
way as the above mentioned windows. For T > 5, the 
precession slows, and the streams settle at a specific loca- 
tion about 30° downstream of (counter-c lockwise from) 
the (Jl, fi) plane. 

Earlier, we noticed that the streams may follow the 
longer path along the dipole magnetic field to the star 
(K02). Now we understand that this is one of the stages 
of the precession around the z— axis. Precession slows 
down because inner regions of the disk come to co- rotation 
with the magnetosphere (see also §5). The precession may 
start again if new matter comes to the inner regions of the 
disk with higher angular momentum. We suggest that if 
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Fig. 8. — Matter flow around the dipole rotator for an inclination angle Q 
isocontours at p = 0.4. 



30°. The background represents the density 



the accretion rate is variable, then this will lead to higher 
angular velocity of precession (at larger M) and smaller 
angular velocity (at lower M). 

Figure 6 shows results of simulations for a larger incli- 
nation angle, 8 = 15°. One can see that matter again 
flows to the star along two streams. However, the streams 
form more rapidly compared to the = 5° case, in 
T « 1.5 rotations. They precess around the z— axis during 
r = 1.5 — 4. Later, they settle into a stationary config- 
uration and impact the star at a location downstream or 
counter-clockwise from the (rj, fi) plane. Figure 7 shows 
four projections at T = 6. Projection (c) shows that the 



streams settled at approximately 30° downstream of the 
(r2, fi) plane, the same as in the = 5° case. 

Analysis of forces along the streams has shown that the 
main lifting force drugging matter up to the streams is the 
pressure gradient force, while the gravity force is the main 
one dragging matter along the stream to the star (see more 
detailed analysis in R03). The magnetic force is negligibly 
small. Similar forces were observed in analysis of funnel 
flows in axisymmetric simulations (R02). However, the in- 
clination of the dipole is favorable for lifting: the lifting 
force can be smaller than in the axisymmetric case. 

The favored locations of the streams on the magneto- 
sphere can be understood as follows. Matter in the inner 
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Fig. 9. — The figure shows different projections for an incUnation angle 8 = 30° and time T — 7. The background 
represents density at the level p = 0.35, and the solid lines are magnetic field lines. In panel (a) the magnetic moment fi 
is directed towards the observer; in panel (b) the orientation is the same as in Figure 8; panel (c) shows the view along 
the Z axis with fi directed downward; and panel (d) shows the side view along the Y axis. 



regions of the disk spirals gradually inward along almost 
circular orbits. When it comes closer to the magneto- 
sphere, it is subject to magnetic braking. The magnetic 
braking is the smallest in the (X, Z) [or, {ft, fi)] plane, 
where the dipole field lines extend to the highest and low- 
est distances above and below the disk. In the {Y, Z) plane, 
the magnetic braking is the largest, because the field lines 
symmetrically cross the disk at this plane. Thus, matter 
flows with low braking through the {X, Z) plane, but is 
slowed down in the (Y, Z) plane. From the other side, the 
poles are closer to the disk matter in the {X, Z) plane, and 
magnetic field lines are directed to the poles, so that only 
small lifting force would be sufficient for the direct accre- 
tion to the pole. In the {Y, Z) plane, matter has much 
longer path to the poles, and larger lifting force should 
be applied for accretion in the (Y, Z) plane. Finally, it 
seems that matter "prefers" to flow in between of {X, Z) 
and [Y, Z) planes. It is important to note that for the 
considered slowly rotating star, the inner part of the disk 
rotates faster than the magnetosphere. Thus the stream 
settles downstream of (counter-clockwise from) the {fl, fi) 
plane. Note, that the location of the stream is different 
for different values of fi* (see R03). 

Figure 6 shows that for T ;> 4, magnetic braking leads 



to a partial disruption of the disk at r < rtr ~ 2—3. In this 
region of the disk, the density becomes lower with time. 
From the other side, matter accumulates near the magne- 
tosphere boundary (r™ « 1 — 1.5) forming a "ring", where 
the density of the disk increased. Matter from the ring 
gradually accretes to the star. This matter distribution 
is similar to that observed in two-dimensional simulations 
of R02. The difference is that in 3D simulations matter 
often forms a spiral structure in the area of the magnetic 
braking r < rhr- 

Figures 6 and 7 show that the inner region of the disk 
becomes asymmetric and the projection (d) shows that in- 
ner parts of the disk are "lifted" above the Z — plane of 
the distant disk. Observations of different cases show that 
this warp is connected with the tendency of the accreting 
matter to co-rotate with the magnetosphere. The inner 
region of the disk is tilted so as to have its normal vector 
in the direction of the magnetic moment fj,. This direction 
of the warp is different from that predicted theoretically 
(see §1). 

Simulations at small inclination angles have shown that 
even a small inclination of the dipole moment leads to non- 
axisymmetric accretion and possible variability of the light 
emitted from the surface of the star. 
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Fig. 10. — Matter flow to an inclined rotator for an inclination angle Q — 45° at the times indicated. The background 
represents the density surface p = 0.35. 



4.2. Medium Inclination Angles: Q = 30° and Q = 45° 

Figure 8 shows the matter flow for an inclination an- 
gle = 30°. We observed that the dense funnel flow 
streams form faster, after about one rotation. At first, 
two prominent streams form and precess as in the case of 
small inclination angles. However, at T = 3 — 4 the mag- 
netic braking starts to disrupt the inner regions of the disk 
(r ^ 2 — 3), and the streams turn into leading spirals. We 
do not bring new matter with Keplerian angular momen- 
tum, so that the density in the inner regions of the disk 
gradually decreases and matter starts to rotate with sub- 
Keplerian velocity. At the same time, a dense ring forms 



around the magnetosphere of the star as in the cases of 
smaller <d. Later, at T ^ 6, the density decreased further 
and a wide gap formed between the inner disk and the 
star. There is still low density matter between the inner 
ring and the inner regions of the disk. This structure is 
similar to that observed in two-dimensional simulations: 
an inner ring forms around the magnetosphere, and there 
is a lower density disk around it. 

Figure 9 shows different projections at T = 6. Projec- 
tion (d) shows that the inner region of the disk is warped. 
The warp arises from the tendency of the accreting matter 
to corotate with magnetosphere; that is, the normal to the 
warped part of the disk tends to align with /i,. 
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Fig. 11. — Different projections of the accretion flow for inclination angle O — 45° and time T = 7. The background 
represents the density surface p = 0.35, and the solid lines are magnetic field lines. In panel (a) the magnetic moment fj, 
is directed towards the observer. Panel (b) shows one interesting projection; panel (c) shows the view along the Z axis 
with /I directed downward, and panel (d) shows the the side view along the F— axis. 



Figure 10 shows the accretion flow at Q — 45°. At 
this relatively large inclination angle, the accreting matter 
forms a complicated structure with initially two streams as 
in other cases. But later, four streams form. Inner regions 
of the disk experience magnetic braking and the disk is sig- 
nificantly disrupted at the chosen density level. However, 
there is always a bridge of matter connecting external re- 
gions of the disk with magnetospheric fiow. This bridge is 
connected with the fact that some matter accretes to the 
poles "directly" as a result of significant inclination of the 
dipole. Such multiple streams may give rise to complicated 
variability of emissions from the stellar surface. 

Two sets of magnetic field lines are shown in Figure 10. 
One set is the same as in Figures 5, 6 and 8. They thread 
the disk at r « 2.8 and are closed inside the region r < 3. 
Another set of field lines start closer to the poles. Some 
of these lines thread the inner regions of the disk due to 
high inclination of the dipole; others thread the disk far 
away, r >> 3. The field lines threading the disk at close 
distances, become involved in the disk rotation and form 
azimuthal magnetic field inside the disk. 

Figure 11 shows different projections of the flow for 
6 = 45° at time T — 7. The projection (a) shows that 
the flow is symmetric if the observer looks along the fi 



axis. Projection (c) shows that the streams associated 
with "direct" accretion to the poles are located at the an- 
gle ~ 20° — 30° downstream (counter- clockwise) from the 
(r2, fi) plane as in all above considered cases with smaller 
6. Projection (d) shows that magnetospheric flow has a 
tendency to co-rotate with magnetosphere as found earlier 
for smaller 0's. The funnel flows rise high above the disk 
plane. 

4.3. Large Inclination Angles: O — 60° and Q = 75° 

Figure 12 shows the evolution of the matter flow for 
8 = 60°. The flow becomes complicated after T > 1. 
From one side, matter accretes directly along field lines to 
the pole which (at such inclination) is close to the plane 
of the disk. From the other side, the matter goes around 
the magnetosphere and accretes through the remote path 
along the dipole field lines. The combination of these two 
flows gives a complicated overall flow. At T = 7, 8, the 
direct polar fiow predominates. However, at T = 9, more 
matter came from the disk and matter flow becomes more 
complicated. 

Two sets of magnetic field lines are shown in Figure 12. 
The number of "polar" field lines threading the inner disk 
is larger than that at 9 = 45° due to the dipole inch- 
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Fig. 12. — Accretion flow to an inclined dipole rotator for an inclination angle Q — 60°. The background represents the 
density surface p — 0.35. 

nation. The magnetic field lines are dragged by the disk closed field lines cross the path of the disk matter, fn such 
forming azimuthal field inside the disk. a configuration, the disk matter experiences strong mag- 
Figure 13 shows four projections of the O ~ 60° flow at netic braking, and it accretes to the poles starting from 
T — 8. One can see that the flow looks very different at dif- T < 1. Some matter, however, tries to continue its "cir- 
ferent projections. Projection (c) shows that most of mat- cular" path around the star. This matter is lifted up and 
ter accretes through the nearby pole of the inclined dipole. goes around the magnetosphere, because the closed mag- 
The location of the streams is slightly downstream of the netosphere is on the way. The magnetic field lines which 
(il, /x) plane. Complicated warping is observed which is initially were closed inside r < 3, moved closer to the star 
connected with magnetospheric accretion. dragged by accreting matter as for other 0's. The "polar" 
Lastly we discuss the high inclination angle flow with magnetic field lines show complicated bending. Some of 
Q = 75°. In this case many "polar" field lines thread the them are dragged by the disk, forming an azimuthal field 
inner region of the disk (see Figure 14 at T = 0) , while the inside the disk. 
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Fig. 13. — Different projections of the accretion flow for an inclination angle 8 = 60° at time T = 8. The background 
represents the density surface p = 0.35; the solid lines are magnetic field lines. In panel (a) the magnetic moment /x is 
directed towards the observer; panel (b) shows one of the interesting projections; panel (c) shows view along the Z axis 
with fjb directed downward; and panel (d) shows the the side view along the K— axis. 



Figure 15 shows projections of the matter fiow at T = 5. 
Projections (a) and (d) show that matter flowing around 
the magnetosphere is strongly lifted above the disk plane. 
These features of the matter flow are typical for a slowly 
rotating star, where matter in the inner regions of the 
disk rotates at a higher angular speed than the star and 
its magnetosphere. Projection (c) shows that the matter 
streams reach the star downstream of the (f2, /x) plane as 
seen in other cases. We conclude that this position of the 
streams at the magnetosphere is a general feature during 
accretion to a slowly rotating star (see also panels (c) in 
Figures 7, 9, 11, 13, and 15). Note that in cases of a faster 
rotating star (closer to equilibrium) or in the case of a very 
fast rotating star ( "propeller" regime) , the location of the 
streams and their shapes are different (see R03). 

4.4. Warping of the Disk 

In all of our simulations we observed that the in- 
ner regions of the disk are disturbed and become non- 
axisymmetric. Some parts of the disk are above and oth- 
ers are below the equatorial plane so that the inner disk 
is warped. In all cases we observed that this warp reflects 
the tendency of matter to co-rotate with magnetosphere 
of the star, so that the normal of the inner warped disk 
has tendency to coincide with the dipole moment jj,. This 



direction of warping does not coincide with any of the pro- 
posed theories. Such a warp is specifically noticeable at the 
Figure 9d, for Q = 30°. This warp has a different nature. 
It represents the matter flow around the magnetosphere 
and is not a part of the disk. One can see that matter is 
warped in the direction of the remote pole. 

A possible reason that the observed warp does not co- 
incide with the theory is that most of theories deal with 
pure diamagnetic disks (e.g., Aly 1980; Lipunov & Shakura 
1980). Lai (1999) considered non-diamagnetic and par- 
tially diamagnetic disk, but his lifting and precessional 
forces mainly rely on the twisting of the magnetic field 
lines, that is on formation of strong component above 
the disk. We did not observe significant twisting of the 
field lines in the inner regions of the disk or in the mag- 
netosphere. In the regions of the disk very close to the 
magnetosphere, <^ r ^ 2 — 3, the disk has tendency 
to corotate with the star. Also, the corona just above the 
disk has a tendency to corotate with the disk (see also 
detailed analysis in R02). For this reason the twisting of 
the magnetic field is small. At larger distances, r > 3, the 
azimuthal component of the magnetic field may be com- 
parable to the poloidal component. However, these field 
lines have tendency to be open. Further, the magnetic 
field in this region is not strong enough to influence the 



17 



T = 




Fig. 14. — Accretion flow to an inclined dipole rotator for an inclination angle O 
density surface p = 0.45. 



75°. The background represents the 



dynamics of the disk. We conclude that the magnetic field 
hues threading the inner regions of the disk are not twisted 
sufficiently to give this type of warping. Further analysis 
is needed for understanding of warping and asymmetries 
observed in simulations. Special simulations with a thin 
diamagnetic disk may be used to check the theory. 

5. ANALYSIS OF THE MAGNETIC BRAKING 

The simulations show that the interaction of the star's 
magnetic field with the inner regions of the disk causes 
a redistribution of density in this part of the disk. The 
three-dimensional pictures presented above show the mag- 
netospheric flow at fixed density levels. To show the spa- 



tial distribution of density, we use two-dimensional slices. 
We also show one-dimensional distribution of density and 
angular momentum density along different axes. 

To compare results obtained at different inclination an- 
gles 8, we take a fixed simulation time. We were able to 
run some cases longer than others, reaching up to 10 — 12 
rotations in some cases. However, in order to incorporate 
the case with the highest inclination, — 75°, we chose 
time T = 5 as the reference time for comparing results at 
different Q. 

Figures 16-18 show two-dimensional slices of the den- 
sity distribution for different inclination angles. Figure 16 
shows the slice of the {X, Z) plane (at Y — Q), which co- 
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Fig. 15. — Different projections of the accretion flow for an inclination angle O = 75° and time T = 5. The background 
represents the density surface p = 0.45; the solid lines are magnetic field lines. In panel (a) the magnetic moment /j, is 
directed towards the observer; panel (b) shows a projection similar to one in Figure 16; panel (c) shows the view along 
the Z axis with jj, directed downward;and panel (d) shows the the side view along the F— axis. 



incides with the plane containing the magnetic moment 
(the ($7, /x) plane). From Figures 16 and 17 it is seen that 
the external regions of the disk, at r ^ 3 — 4, are not ap- 
preciably disturbed and are approximately the same for 
different inclination angles 0. Inside this region, however, 
the scales at which the magnetosphere influences the mat- 
ter flow are different. Figure 16 shows that for G = 5°, 
matter starts to go into funnel flows around the magne- 
tosphere at approximately r w 1.5. For 8 — 15°, this 
occurs at r « 1.8; at 9 = 30°, at r « 1.9. At 6 45°, 
60°, and 75°, the region is even larger, r k, 2.2. Fig- 
ure 17 shows the density distribution in the (y, Z) plane 
(at X — Q). This cross-section also shows that at larger 
8, the magnetospheres are larger. Figure 18 shows the 
(X, Y) cross-sections (Z = 0), which go through the plane 
of the disk. Figure 18 shows that a relatively high density 
ring (white color) forms around the low-density magneto- 
spheres at small inclination angles, 9 = 5°, 15° and 30°, 
but it disappears at larger inclination angles. At very large 
inclination angles, Q = 60° and 75°, the direct accretion 
to the poles is observed (white color in the polar regions). 
This cross-section also shows that the inner region with 
radii i? 3 — 4, has a lower density than the rest of the 
disk. 



We also did one- dimensional analysis of the density and 
angular momentum distribution along the X and Y axes. 
For clarity of the presentation, we divided the inclination 
angles into two groups: smaller angles, 6 = 0°, 5°, 15°, 
and 30°, and larger angles: 6 = 45°, 60°, 75°. The case 
9 = 30°, was added to the second group for convenience of 
comparisons. Figure 19a shows density distribution along 
the X axis for smaller inclination angles. One can see 
that the disk is disturbed by magnetic braking for r < A. 
One can see that at the larger inclination angle the den- 
sity is smaller. The density is very small for X <^\^ where 
the magnetic field is very strong. The strong magnetic 
field limits the matter flow. At larger X, the density in- 
creases forming a peak of density. This is the ring men- 
tioned earlier, where matter accumulates before accreting 
to the poles. This feature was the typical one in all two- 
dimensional simulations (see R02). It is also clearly ob- 
served in the simulations at small inclination angles. Fig- 
ure 19b shows similar distributions along the Y axis, which 
are similar to the distributions along the X axis. 

At larger inclination angles, 8 = 45°, 60° and 75°, the 
density distribution is different (see Figures 19c, d). In the 
X direction, for 8 = 45°, matter comes closer to the stel- 
lar surface. For 9 = 60° and 8 = 75°, the density strongly 
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Fig. 16. — The (X, Z) plane cross-section (in which magnetic moment is located) showing the density distribution at 
T = 5 for different inclination angles. 



increases towards the star, which means that matter ac- 
cretes directly to the star. The magnetoshere is not an 
obstacle for the flow at high inclination angles 8. In the 
Y direction (see Figure 19d), the matter is stopped by the 
magnetospheric at a radius r « 1 for all inclination angles. 

Magnetic braking leads to the deviation of the angular 
velocity of the matter from the Keplerian value. Figures 
20a, b show that at radii "r ^ 3, the angular velocity uj devi- 
ates from Keplerian value starting from r <; 3. It gradually 
become smaller than Keplerian and at r <J 1.3 it sharply 
decreases to the angular velocity of the star SI, — 0.19. At 
large inclination angles (see Figures 20c, d) the deviation 
from Keplerian rotation is much stronger in both the X 



and Y directions. For 8 = 60° and 75°, the angular ve- 
locity inside the region r <; 3 is w ~ 0.2 — 0.3; that is, the 
inner regions of the disk almost co-rotate with the star. 

6. MATTER AND ANGULAR MOMENTUM FLUXES 

Here, we analyze the mass accretion rate to the surface 
of the star for different inclination angles. Also, we ana- 
lyze the angular momentum flux to the star and discuss 
its influence on the angular momentum of the star. 

6.1. Matter Flux 
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Fig. 17. — The {Y, Z) plane cross-section showing the density distribution for T = 5 at different incUnation angles. 



We calculated the matter flux to the surface of the star, 



M 



dS ■ pv , 



(5) 



where the surface integral is taken just outside of the star 
with dS the outward pointing surface area element. 

Figure 21 (left panel) shows the matter fluxes at differ- 
ent inclination angles versus time. Initially the accretion 
rate is zero, but it increases rapidly as matter reaches the 
surface of the star. It takes less time for larger inclination 
angles. After the initial rise, the accretion rate settles at 
a value which changes only slowly, forming the "plateau" . 
At = 30° and 45°, the plateau values of M are approxi- 
mately the same, M « 0.19. At 6 = 60°, M « 0.2 - 0.24, 



and at e = 75° it is even larger, M « 0.26. At the 
smallest inclination angles, = 5° and = 15°, the M 
increases gradually and reaches the values M « 0.19 — 0.2, 
but it is not clear, whether the plateau was formed or 
not. Note, that at Q = 0°,the accretion rate is smaller, 
M w 0.13, but not significantly smaller than that for in- 
clined rotators. These values of M are in accord with two- 
dimensional axisymmetric simulations done with a similar 
grid, M Ki 0.14 — 0.16. Note, that simulations in the hy- 
drodynamic case {B — 0) lead to a much smaller accretion 
rate M « 0.02. This accretion is due to the small nu- 
merical viscosity as discussed in the Section 2. Thus, we 
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Fig. 18. — The {X, Y) plane cross-section showing the density distribution bXT — b for different inclination angles. 



conclude that the accretion rate is larger at larger 8 be- 
cause the magnetic braking is larger at larger Q. 

6.2. Torque on the Star 

Incoming matter and the magnetic field carry in or re- 
move angular momentum from the star. The angular mo- 
mentum of this matter is much less than that of the star 
during the time span of the simulations. Thus, we fix an- 
gle O and direction and value of Jl. We calculate torques 
relative to different axes which carry information about 
possible long-term evolution of the star. 

The total angular momentum flux in the reference frame 
rotating with the star {X, Y, Z) consists of the part car- 
ried by the matter which gives the matter torque, Nm, 



and that carried by the field which gives the field torque 
N/, 



J dS-v p (r X u)+^ J dS-B (r x B), 



N = N,„+N/ = - 

where the integration is over a spherical surface just out- 
side of the star with dS the outward pointing surface area 
element, and where u is the velocity in the inertial frame 
and V is the velocity in the frame rotating with the star. 
For an inclined rotator, N is not in general along a partic- 
ular axis. Denoting L the angular momentum of the star, 
we have 

dL d*L ^ ^ 

rJxL = N (7) 



dt 



dt 
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Fig. 19. Density distribution along the x (left panels) and y (right panels) axes for inclination angles G = 0°, 5°, 15°, 30° 
(top panels) and = 30°, 45°, 60°, &75° (bottom panels). 



where d/dt is with respect to an inertial reference frame 
and d*/dt is with respect to the rotating {X,Y,Z) coor- 
dinate system. For the considered case, the angular mo- 
mentum of the star and that of the incoming matter in the 
disk have the same direction. Therefore, the direction of 
L averaged over the star's rotation period remains in the 
z direction of the inertial or laboratory reference frame. 

We observed that in the vicinity of the star the angular 
momentum flux is carried predominantly by the magnetic 
field. The flux carried by the matter is very small, about 
1% of the total. Far from the magnetosphere, the angu- 
lar momentum flux is carried mainly by the matter. With 
decreasing distance, the matter transfers its angular mo- 
mentum to the magnetic field. A similar behavior was 
observed in axisymmetric simulations (R02). 

The torque about the Z axis, N^, is responsible for 
spinning-up or spinning-down of the star. Positive cor- 
responds to spinning-up of the star. Figure 21 (right panel) 
shows the flux Nz for different inclination angles. In the 
present work the star rotates slowly so that the incoming 
matter spins-up the star. The variation of is analogous 
to the variation of M (Figure 21. left panel). One can 
sec, that initially, the magnetic braking was large for all 
inclination angles, but later it decreased and became rel- 
atively constant. At small inclination angles, Q = 0°,5° 
and 15°, the torque went to an approximately constant 
value rapidly, T < 0.5. The torque is almost constant for 
= 0°,5° , and decreases only slowly at O = 15°. At 
large 0, the decreases gradually. If new matter is ac- 
creted to the inner regions of the disk, then A^^ will stay at 



some non-zero level, as we observed in our two-dimensional 
simulations (R02). Additional simulations with steady ac- 
cretion from the disk are needed to resolve the difference 
in angular momentum transport at different 6. 

A positive torque would, in the absence of the $7 x L 
term, act to shift the star's angular momentum L in the 
+X direction, that is towards the direction of fi. How- 
ever, the r2 X L is expected to be important for an actual, 
oblate star where the moments of inertia are Iz > Ix = ly 
with {Iz — Ix)/Ix ^ 1. For this reason a torque Nx acts 
to give Qy — Nx/ (Ix^z^prec) which corresponds to a very 
slow precession of the $7 about L which is fixed in space, 
where flprec = {Iz — Ix)^z/Ix is the free precession an- 
gular frequency. Figure 22 (left panel) shows the torque 
Nx for different inclination angles 0. One can see, that at 
= 0, Nx = 0, which confirms that our three-dimensional 
code supports axisymmetry with high accuracy. This also 
confirms that the origin of this momentum is connected 
with interaction of the inclined dipole with the disk. Even 
at small inclination angle, = 5°, the torque A^^ is non- 
zero. It varies sinusoidally and changes sign from negative 
to positive and back, and settles to a positive value. For 
= 30° , variations of larger amplitude were observed. At 
inclination angles, = 45° and 60°, Nx is initially neg- 
ative, but becomes positive later. A negative torque at 
= 75° gradually increases to larger values. 

A positive torque Ny would, in the absence of the x L 
term, act to shift the star's angular momentum L in the 
+Y direction, that is towards the direction of x /it. How- 
ever, as mentioned the 12 x L is expected to be important 
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Fig. 20. — Angular velocity distribution along the X axis (left panels) and Y axis (right panels) for inclination angles 
e = 0°, 5°, 15°, 30° (top panels) and 6 = 30°, 45°, 60°, & 75° (bottom panels). 
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Fig. 21. — The left-hand panel shows the matter flux to the star M versus time T for different inclination angles O. The 
right-hand panel shows torque versus time at different G. 



for an actual, oblate star. For this reason a torque iVj^ 
acts to give 0,^ = —Ny/{Ix^zQprec) which corresponds to 
a very slow precession of the CI about L which is fixed in 
space. Figure 22 (right panel) shows the torque A^j^ for 
different inclination angles. We observed that for = 0, 
Ny = 0, as it should be. At the small inclination angle 



e = 



N„ varies around zero value, which is connected 



with the formation of two streams which precesses around 
the magnetosphere of the star. At G = 15°, the sign of Ny 
is mainly negative, though variations of the flux arc also 
observed. At G = 30° and G = 60°, the initial variation of 
the flux at T 2 has the same value, Ny —0.18. How- 



ever, later both fluxes decreased to smaller values. Note 
that at G = 30°, the Ny changes sign to positive values. 
For G = 75°, the torque Ny goes to an approximately 
constant value w —0.07. 

We emphasize that for more rapidly rotating stars, the 
torques may behave quite differently from the case of the 
slowly rotating star considered here. 

7. CONCLUSIONS AND OBSERVATIONAL CONSEQUENCES 

We performed full three-dimensional ideal MHD sim- 
ulations of disk accretion to a slowly rotating star with 
an inclined dipole magnetic field for different inclination 
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Fig. 22. — Torque on the star along the X axis (left panel) and Y axis (right panel) at different inclination angles Q. 



angles G. Below, we summarize the main results of the 
simulations (§7.1) and discuss possible observational con- 
sequences (§7.2). 

7.1. Main Results of Numerical Simulations 

In the following we summarize the main results and con- 
clusions of our simulations. 

1. We observed that the accretion flow near the star be- 
comes non-axisymmetric at very small inclination angles, 

e - 2° - 5°. 

2. The shape of magnctosphoric flow depends on the 
density. At a sufficiently low density, most of the volume 
of the magnetosphere is filled by matter, excluding small 
"windows" of even lower density. At a higher density, the 
matter is in two or several "accretion streams." 

3. For relatively small inclination angles, G <; 30°, 
the densest matter accretes in two streams. The streams 
typically follow a path to the closest magnetic pole of the 
star with displacement of 20° — 30° downstream (counter- 
clockwise) from the ($7, /j,) plane on the star's surface. 

The accretion streams may precess about the star's rota- 
tion axis. The precession appears to be caused by incom- 
ing matter with higher than average angular momentum; 
this matter causes the streams to rotate faster than the 
star. In periods of steady accretion, the streams do not 
precess. Precession is observed more often in cases with 
small G. 

For intermediate inclination angles, 30° < G < 60°, 
the magnetospheric flow consists of two or more streams. 
Several streams (typically four) is a common feature for 
G = 45°, while in other cases (G = 30°, G = 60°) sev- 
eral streams appear in periods of enhanced accretion. The 
number of streams may depend on the density level. 

For large inclination angles, G ^ 60°, matter typically 
accretes in two streams. However, the streams have a dif- 
ferent shape compared with those at small G and they 
come to the star near the equatorial plane. 

4. The density distribution is close to that observed 
in two-dimensional simulations (R02). Namely, magnetic 
braking leads to depression of density in the region ;^ 
r and to the formation of a dense ring near the mag- 
netospheric boundary, r ~ r^. Formation of a dense ring 
is typical for G <^ 30°. At larger G, the ring does not form, 
and for G ^ 60°, matter accretes to the star through equa- 
torial streams (see Figures 16-18). 



Magnetic braking leads to a significant departure of the 
angular velocity of the disk from Kcplcrian in the region 
r ^ 2 - 3 For G = 60° and G = 75°, the region of the disk 
out to r w 2.5 almost co-rotates with the star (see Figure 
20c,d). 

5. The inner regions of the disk are often warped. This 
warp forms because matter leaving the disk, start to move 
around magnetosphere, thus forming a flow with axis close 
to that of n- The region between the undisturbed disk 
and the matter which co-rotates with magnetosphere has 
a warped shape. However, it has a different nature com- 
pared to that predicted by theory of warped disks (see 
§4.4). 

6. The angular momentum transport to the star was 

investigated. For the considered small angular velocity of 
the star, the angular momentum flux relative to the Z 
axis gives a positive torque which acts to spin-up of 
the star for all G. Angular momentum is transported to 
the surface of the star by the magnetic field (as observed 
in our two-dimensional simulations, R02) The matter car- 
ries about 1% of the total angular momentum fiux near 
the surface of the star. 

The torques and Ny were calculated in the reference 
frame corotating with the star. They give information 
about the long term evolution of the star's angular mo- 
mentum L, and of and /ix. For typical accretion rates 
(~ 10~^— lO~^M0/yr for a protostellar system), the angu- 
lar momentum of the star changes by a neglibiblc amount 
during a rotating period of the inner region of the disk. 
This justifies our having G fixed in the simulations. 

7.2. Observational Consequences 

The described simulations are important for under- 
standing the physics of accreting magnetized stars, for ex- 
ample, CTTSs, cataclysmic variables, and X-ray pulsars. 
The simulations provide a basis for understanding the pho- 
tometric and spectral variability the stars. Radiation in 
different spectral bands may be associated with different 
regions, the stellar surface, hot/cold spots on the stellar 
surface, and the magnetospheric plasma. Knowledge of the 
magnetospheric flow is essential to interpreting observed 
light curves, spectra, amd temporal variations in the spec- 
tra. The presented simulations may allow the estimation 
of G in different systems, and may help to establish the 
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main source of the radiation in different wavebands. Some 
of predictions from this work are: 

1. Non-axisymmetric accretion and the formation of ac- 
cretion streams (which give hot spots on the star's surface) 
is expected in almost all stars, because the limiting angle 
Q at which non-axisymmetry appears is small {Q ~ 2°). 

2. At a fixed G, the geometry of the magnetospheric 
flow is different at different density levels. At a sufHciently 
low density, all space is occupied by matter at this or larger 
densities. At a larger density level, "windows" appear in 
the accretion flow, and at even larger densities, matter ac- 
cretes along narrow accretion streams. The variability of 
different spectral lines and of different wavebands of the 
continuum radiation is expected to depend on the den- 
sity, temperature and velocity distributions along the line 
of sight to the star. Variation of the accretion rate may 
also change the pattern of variability, because at the same 
O the number of streams may temporary change. Spiral 
structure may appear in the inner regions of the disk. 

3. If the light curve has features with variability equal 
to approximately ahalf of the main period, this may be a 
sign of two rotating hot/cold spots, or, a sign of occupa- 
tion of stellar light by two magnetospheric streams. If the 
light curve suggest occultations, then the inclination angle 
predicted by simulations is relatively small, O < 30°. Rel- 
atively short periods (few rotations) of the multi-stream 
features may appear, but most of the time the two-stream 
features will dominate. Different observations point to 
such a possibility (e.g., MuzeroUe, Hartmann & Calvet 
1998, Petrov et al. 2001). 

4. If the light curve indicates precession of the accre- 
tion streams, then the simulations point to the likelihood 



that the inclination angle is small, O <; 10°. The variabil- 
ity pattcaii for two precessing streams is characterized by 
an initial period smaller than the half-period P/2, and a 
subsequent regular increase of period to P/2. This type of 
variability may be rare because at small inclination angles, 
the expected amplitude of variability is small. 

5. If the light curve shows a quasi- period much shorter 
than the period of rotation of the star (say, ~ l/AP), 
then it will be a sign of multiple streams which are typical 
for intermediate inclination angles, 30° < 9 < 60°. Re- 
cently MuzeroUe, Calvet, & Hartmann (2001), suggested 
that some observations of CTTSs may be explained if the 
funnel flows have a multiple stream geometry. 

6. At large inclination angles, 8 ^ 60°, matter also 
accretes in two streams, but they arc located close to the 
equatorial plane. In such a case variability may be con- 
nected with occultation of stellar light by the warped disk 
(as proposed by Bouvier et al. 1999, 2003). The light 
curves are expected to be different from those where oc- 
cultation by magnetospheric streams dominate. 
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